#ifndef MATRIX_HPP
#define MATRIX_HPP

#include <vector>
#include <stdexcept>
#include <iostream>

template<typename T>
class Matrix {
public:
    // 构造函数
    Matrix() : rows(0), cols(0) {}

    Matrix(size_t r, size_t c) : rows(r), cols(c), data(r * c, T()) {}

    // 拷贝构造函数
    Matrix(const Matrix& other) : rows(other.rows), cols(other.cols), data(other.data) {}

    // 静态工厂方法
    static Matrix Zero(size_t r, size_t c) {
        return Matrix(r, c);
    }

    static Matrix Identity(size_t size) {
        Matrix m(size, size);
        for (size_t i = 0; i < size; ++i) {
            m(i, i) = T(1);
        }
        return m;
    }

    // 访问元素
    T& operator()(size_t i, size_t j) {
        if (i >= rows || j >= cols) {
            throw std::out_of_range("Matrix index out of range");
        }
        return data[i * cols + j];
    }

    const T& operator()(size_t i, size_t j) const {
        if (i >= rows || j >= cols) {
            throw std::out_of_range("Matrix index out of range");
        }
        return data[i * cols + j];
    }

    // 矩阵运算
    Matrix operator*(const Matrix& other) const {
        if (cols != other.rows) {
            throw std::invalid_argument("Matrix dimensions don't match for multiplication");
        }

        Matrix result(rows, other.cols);
        for (size_t i = 0; i < rows; ++i) {
            for (size_t j = 0; j < other.cols; ++j) {
                T sum = T();
                for (size_t k = 0; k < cols; ++k) {
                    sum += (*this)(i, k) * other(k, j);
                }
                result(i, j) = sum;
            }
        }
        return result;
    }

    Matrix operator+(const Matrix& other) const {
        if (rows != other.rows || cols != other.cols) {
            throw std::invalid_argument("Matrix dimensions don't match for addition");
        }

        Matrix result(rows, cols);
        for (size_t i = 0; i < rows; ++i) {
            for (size_t j = 0; j < cols; ++j) {
                result(i, j) = (*this)(i, j) + other(i, j);
            }
        }
        return result;
    }

    Matrix operator-(const Matrix& other) const {
        if (rows != other.rows || cols != other.cols) {
            throw std::invalid_argument("Matrix dimensions don't match for subtraction");
        }

        Matrix result(rows, cols);
        for (size_t i = 0; i < rows; ++i) {
            for (size_t j = 0; j < cols; ++j) {
                result(i, j) = (*this)(i, j) - other(i, j);
            }
        }
        return result;
    }

    // 转置
    Matrix transpose() const {
        Matrix result(cols, rows);
        for (size_t i = 0; i < rows; ++i) {
            for (size_t j = 0; j < cols; ++j) {
                result(j, i) = (*this)(i, j);
            }
        }
        return result;
    }

    // 求逆（使用高斯-约旦消元法）
    Matrix inverse() const {
        if (rows != cols) {
            throw std::invalid_argument("Matrix must be square for inverse");
        }

        // 创建增广矩阵 [A|I]
        Matrix augmented(rows, 2 * cols);
        for (size_t i = 0; i < rows; ++i) {
            for (size_t j = 0; j < cols; ++j) {
                augmented(i, j) = (*this)(i, j);
                augmented(i, j + cols) = (i == j) ? T(1) : T(0);
            }
        }

        // 高斯-约旦消元
        for (size_t i = 0; i < rows; ++i) {
            // 找主元
            T pivot = augmented(i, i);
            if (pivot == T(0)) {
                throw std::runtime_error("Matrix is singular");
            }

            // 归一化当前行
            for (size_t j = i; j < 2 * cols; ++j) {
                augmented(i, j) /= pivot;
            }

            // 消元
            for (size_t k = 0; k < rows; ++k) {
                if (k != i) {
                    T factor = augmented(k, i);
                    for (size_t j = i; j < 2 * cols; ++j) {
                        augmented(k, j) -= factor * augmented(i, j);
                    }
                }
            }
        }

        // 提取逆矩阵部分
        Matrix result(rows, cols);
        for (size_t i = 0; i < rows; ++i) {
            for (size_t j = 0; j < cols; ++j) {
                result(i, j) = augmented(i, j + cols);
            }
        }
        return result;
    }

    // 获取维度
    size_t get_rows() const { return rows; }
    size_t get_cols() const { return cols; }

    // Scalar multiplication
    Matrix operator*(T scalar) const {
        Matrix result(rows, cols);
        for (size_t i = 0; i < data.size(); ++i) {
            result.data[i] = data[i] * scalar;
        }
        return result;
    }

    friend Matrix operator*(T scalar, const Matrix& mat) {
        return mat * scalar;
    }

private:
    size_t rows;
    size_t cols;
    std::vector<T> data;
};

#endif 